Quantum algorithms for solving linear differential equations 



Dominic W. Berry 1 

'institute for Quantum Computing, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada 

Linear differential equations are ubiquitous in science and engineering. Quantum computers can 
simulate quantum systems, which are described by homogeneous linear differential equations that 
produce only oscillating terms. Here we extend quantum simulation algorithms to general inho- 
mogeneous linear differential equations, which can include exponential terms as well as oscillating 
terms in their solution. As with other algorithms of this type, the solution is encoded in amplitudes 
of the quantum state. The algorithm does not give the explicit solution, but it is possible to extract 
global features of the solution. 

I. INTRODUCTION 

Differential equations are used for an enormous variety of applications, including industrial design and weather 
prediction. In fact, many of the main applications of supercomputers are in the form of large systems of differential 
equations [l|. Therefore quantum algorithms for solving differential equations would be extraordinarily valuable. A 
quantum algorithm for differential equations was proposed in Ref. [2j, but that algorithm had very poor scaling in 
the time. The complexity of the simulation scaled exponentially in the number of time-steps over which to perform 
the simulation. 

The algorithm in Ref. @| may have been overly ambitious, because it aimed to solve nonlinear differential equations. 
A more natural application for quantum computers is linear differential equations. This is because quantum mechanics 
is described by linear differential equations. We find that, when we restrict to linear differential equations, it is possible 
to obtain an algorithm that is far more efficient than that proposed in Ref. [2j. 

We consider first-order linear differential equations. Using standard techniques, any linear differential equation 
with higher-order derivatives can be converted to a first-order linear differential equation with larger dimension. A 
first-order ordinary differential equation may be written as 

x(t) = A(t)x{t) + b{t), (1) 

where x and b are iV^-component vectors, and A is an N x x N x matrix. Classically, the complexity of solving the 
differential equation must be at least linear in N x . The goal of the quantum algorithm is to solve the differential 
equation in time O (poly log N x ). 

Quantum mechanics is described by differential equations of this form, except they are homogeneous (b(t) = 0), and 
A(t) — iH(t), where H(t) is Hcrmitian. This means that the solutions in quantum mechanics only include oscillating 
terms, whereas more general differential equations have solutions that may grow or decay exponentially. Quantum 
algorithms for simulating quantum mechanical systems have been extensively studied [3|— 19|] . 

Classical physics is described by more general differential equations. Large systems of ordinary differential equations 
are produced by discretisation of partial differential equations. Many equations in physics are linear partial differential 
equations, where the time derivative depends linearly on spatial derivatives and the value of a quantity at some point 
in physical space. Examples include Stokes equations (for creeping fluid flow), the heat equation, and Maxwell's 
equations. Discretisation of the partial differential equation on a mesh of points results in an ordinary differential 
equation with a very large value of N x . 

In the case where A and b are time independent, then one can find the equilibrium solution of the differential 
equation by solving 

Ax = -b. (2) 

A quantum algorithm for this problem was given by Harrow, Hassadim and Lloyd [Io| , with runtime that is polynomial 
in log(A^) and the condition number of A. Ambainis has reported development of an improved algorithm [111 ], though 
this algorithm has not yet been released. We consider the more difficult case of solving the time evolution under linear 
differential equations, rather than just the equilibrium solution. We find that this case can also be solved using a 
modification of the method of Harrow, Hassadim and Lloyd. 

II. TROTTER FORMULA APPROACH 

Before explaining that approach, we first describe an approach using Trotter formulae, and the drawback to that 
approach. This will not be described rigorously, because it is not our main proposal for solving differential equations. 
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The homogeneous case, where b = 0, is analogous to Hamiltonian evolution. If A is antiHermitian, then we can 
take A = iH, where H is a Hermitian Hamiltonian. Evolution under this Hamiltonian can be solved by methods 
considered in previous work [(| Another case that can be considered is where A is Hermitian. In this case, the 
eigenvalues of A are real, and A can be diagonalised in the form A = V DV~ X , where D is a real diagonal matrix and 
V is unitary. The formal solution is then, for A independent of time, x(t) = Ve D ^ to ^ V~ 1 x(to). 

The differential equation can be solved using a similar method to that used in Ref. [lflj. The value of x is encoded 
in a quantum state as 

\x)=M x J2x [j] \3), (3) 

where \ j) are computational basis states of the quantum computer, x^' are the components of the vector x, and Af x 
is a normalisation constant. The state can be written in a basis corresponding to the eigenvectors of A: 

3 

Using methods for Hamiltonian simulation, iA can be simulated. By using phase estimation, if the state is an 
eigenstate then the eigenvalue Xj can be determined. Given maximum eigenvalue A max , we would change the 
amplitude by a factor of e C*-*o)(A i -A m «) < See Ref JJ^ 

for the method of changing the amplitude. If this is done 

coherently, then the final state will encode x(t). 

For more general differential equations, A will be neither Hermitian nor antiHermitian. In this case, one can break 
A up into Hermitian (Ah) and antiHermitian (A a u) components. The evolution under each of these components can 
be simulated individually, and the overall evolution simulated by combining these evolutions via the Trotter formula. 
The drawback to this approach is that it appears to give a complexity that increases exponentially with the time 
interval At = t — to (though the complexity is still greatly improved over Ref. Q). 

If A were just Hermitian, then the eigenvector (or eigenspace) corresponding to the largest eigenvalue would not 
decay, and the system would end up in that state. Therefore the amplitude would not drop below the amplitude on 
the eigenspace corresponding to the largest eigenvalue. That is not the case when A is a more general matrix, because 
usually the maximum real part of an eigenvalue of A will be strictly less than the maximum eigenvalue of Ah- The 
amplitude must therefore decay exponentially, because we must use the maximum eigenvalue of Ah in simulating 
evolution under Ah- 

The result of this is that the complexity of the simulation will scale exponentially in the time that the differential 
equation needs to be simulated over, At. The scaling will be considerably improved over that in Ref. 0, but it 
is desirable to obtain scaling that is polynomial in At. Another drawback is that this approach does not enable 
simulation of inhomogeneous differential equations. 



III. LINEAR SYSTEMS APPROACH 

To avoid this problem we propose an approach based on the algorithm for solving linear systems from Ref. [To| . 
The trick is to encode the solution of the differential equation at different times using the one state. That is, we wish 
to obtain the final state proportional to 

M:=$>i>l*i>- (5) 

3=0 

The number Nt is the number of time steps, t 3 - is the time to + jh, where h is the time interval in the discretisation 
of the differential equation, Xj is the approximation of the value of x at time tj, and At is the total time interval 
over which the differential equation is to be solved. We use the subscript j to index the vectors, and superscript for 
components of these vectors. 

Once this state has been created, the state encoding the solution at the final time to + At can be approximated by 
measuring the register encoding the time and getting that time. Just using this method, the probability of obtaining 
the final time is small (l/(N t + 1)). To obtain a significant probability of success, one can add times beyond to + At 
where x is constant. We take x to be constant for to + At to to + 2 At, so Nt = 2At/h. Then any measurement 
result for the time in this interval will give the state corresponding to the solution. By this method, the probability 
of success can be boosted significantly, without changing the scaling for N t . 
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To numerically solve differential equations, the simplest method is the Euler method, which discretises the differ- 
ential equation as 

Xj+1 ~ Xj =A(t j )x j +b(t j ). (6) 

For times after to + At, we set Xj+i = Xj to ensure that x is constant. The Euler method yields an error that scales 
as 0{h 2 ) for a single time step. Therefore, we expect that the error in the total simulation is 0(N t h?) = 0(At 2 /N t ). 
To achieve error bounded by e, we can take N t — 0(At 2 /e). To show these scalings rigorously requires additional 
constraints on the problem. 

In particular, to rigorously bound the error it is necessary that the eigenvalues of A(tj) have no positive real part. 
Otherwise the error can grow exponentially. In cases where A(tj) does have an eigenvalue with positive real part, one 
can simply subtract a multiple of the identity, and rescale the solution. Note that e is the error in the solution of the 
differential equation, and is distinct from error in the solution of linear systems. 

More generally, linear multistep methods have the form [l2l . [l3| 

k k 

aexj+e = hJ2/3 e [A(t J+ e)x J+e + b(t j+i )]. (7) 

£=0 i=0 

Multistep methods can be chosen such that the error is of higher order in h, but there is the problem that the method 
may not be stable. That is, even if the exact solution of the differential equation is bounded, the solution of the 
difference equation may be unbounded. 

To examine the stability, one defines the generating polynomials 

k k 

p(c) = £ a ^'' <KO = £&c j '. (8) 

3=0 j=0 

The stability can be examined via the roots of the equation 

p(0 - MC) = o. (9) 

One defines the set S by 

s ■= lu g c- a11 r00ts of ® satisf y I0(m)I < 1 \ ( 10 -) 

V ' multiple roots satisfy | I < 1 J 

S is called the stability domain or stability region of the multistep method. In addition, if the roots of <r(£) all satisfy 
ICI < 1, and repeated roots satisfy |£| < 1, then the method is said to be stable at infinity. 

A linear multistep method is said to be order p if it introduces local errors 0(h p+1 ). This means that, if it is 
applied with exact starting values to the problem x = t q (0 < q < p), it integrates the problem without error. A 
linear multistep method has order p if and only if [l2j 

p{e h ) - ha{e h ) = 0{h p+1 ). (11) 

A useful property of linear multistep methods is for them to be >l-stable [HI, [l4[ • 

Definition 1. A linear multistep method is called A-stable if S D C~ , i.e., if 

ReA < =>■ numerical solution for x = Ax is bounded. (12) 

This definition means that, if the solution of the differential equation is bounded, then the approximation given 
by the multistep method is bounded as well. For a scalar differential equation, the multistep method is bounded 
whenever A is in the left half of the complex plane. The Euler method is ^4-stable, but it is not possible to construct 
arbitrary order A-stable multistep methods. The second Dahlquist barrier is that an A-stable multistep method must 
be of order p < 2 HH3. As we wish to consider hi ghe r-order multistep methods, we relax the condition and require 
that the linear multistep method is A(a)-stable (l3l Il5|. 

Definition 2. A linear multistep method is A(a)-stable, < a < n/2, if 



S D S a = {>;|arg(-/i)| < a,p f 0}. 



(13) 
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This definition means that, in the case of a scalar differential equation, the multistep method is bounded whenever 
A is within a wedge in the left half of the complex plane. For a vector differential equation, the eigenvalues of A should 
be within this wedg e. It is known that, for any a < tt/2 and k G N, there is an j4(a)-stable linear fc-step method of 
order p = k QJ, [iff. 

The error in the total solution of the differential equation will be 0(N t (At) p+1 ). In order to obtain a rigorous 
result, we specialise to the case that A and b are independent of time. The relevant bound is given in Theorem 7.6 in 
Chapter V of Ref. [lj]. 

Theorem 3. Suppose a linear multistep method is of order p, A(a)-stable and stable at infinity. If the matrix A is 
diagonalisable (i.e. there exists a matrix V such that V~ 1 AV = D = diag(Ai, . . . , X n )) with eigenvalues satisfying 

|arg(-A 2 )|<« fori = l,...,N x , (14) 

then there exists a constant M (depending only on the method) such that for all h > the global error satisfies 

\\x(t m ) - x m \\ < Mk v ( max - Xj \\ + [ ™ \\x^ +1 \0\\^) , (15) 

\0<j<k J tQ J 

where Ky — ||V|| • is the condition number ofV. 

Here the superscript with round brackets denotes repeated derivative. We can use this result to show a lemma on 
the scaling of the error. 

Lemma 4. Suppose a linear multistep method is of order p, A(a) -stable and stable at infinity. If the matrix A is 
diagonalisable (i.e. there exists a matrix V such that V~ 1 AV = D = diag(Ai, . . . , X n )) with eigenvalues satisfying 

|arg(-Ai)|<a for i = 1, . . . , N x , (16) 

and b is constant, then the global error satisfies 

\\x(t m ) - x m \\ = O (4(lknit|| + \\b\\/\\A\\) M^IIAII) 2 + m(h\\A\\r +1 ]) , (17) 

where Ky = \\V\\ ■ ||V^ _1 || is the condition number of V . 

Proof. The linear multistep method requires a starting method to obtain the values of Xj for < j < k. The term 
maxo<j<fc || x(tj) — Xj\\ arises from the inaccuracy in this procedure. One can simply use the Euler method, in which 
case the error is 0(h 2 ). It is also possible to use higher-order starting methods, but there is not a convenient rigorous 
result that can be used. To determine the error in the Euler method, one can simply use Theorem |3 Because k = 1 
and p = 1 for the Euler method, and for the initial point there is zero error (x(to) — %o), Theorem [3] gives 

\\x(tj) - xA\ < M E n v h j 3 |M 2 )(OR, (18) 



where Me is the constant for the Euler method. We consider this expression for < j < k, and obtain 

max WxOtA-xA <M E n v h 2 (k-l) max ||a; (2) (£)ll- (19) 

0<]<k ^ e [t 0: t + {k-l)h] 

In using these results, it is necessary to place upper bounds on the values of ||a;' p+1 - ) ('?)ll an( i (£)|| . In general 
these will depend on the value of b(t), and its time-dependence. It is well-behaved if b is a constant, in which case the 
exact solution is 

x(t) = e^-^ix^t + A-H) - A- X b. (20) 

Then 

x^(t) = e A ^(A l x init + A l ^b), (21) 

so 

\\xW(t)\\ - WVe^-^V-^A'x-^ + A^W 

<n v {\\A\\%x iriit \\ + \\A\\^\\b\\) (22) 
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In the first line we have used the diagonalisation of A, and in the second line we have used the condition that 
| arg(— Aj)| < a. Using Theorem [31 the error is bounded as 



ll^(im) — %m\\ < Mkv max ||x(£i) — xA\ + h p 

\0<j<k Jt 



,(p+l) 



(OK 



< Mk v [M E K 2 v h 2 k(\\x- m - lt \\ + ||6||/P||)P|| 2 + hP(t m - *o)«v(||a!init|| + ||6||/P||)||A||P +1 ] 
= O (4(^init|| + \\b\\/\\A\\) M^PII) 2 + ro(ft||A||)*+ 1 ]) . 



(23) 

□ 



This result means that, disregarding the dependence on many of the quantities, and omitting the error due to the 
starting method, the error scales as 0((||^4||ft.) p || A\\ At) for total time At. To achieve error bounded by e, we then use 



N t = 



(\\A\\At) 



l+i/p 



That is, the number of time steps required is close to linear in the time. 

Given a linear multistep method, it is straightforward to encode this method as a linear system 

Ax = b. 



(24) 



(25) 



Here x is the vector of blocks Xj , b is a vector of the blocks b, and A is a matrix describing the initial condition and 
discretiscd differential equation. As an example, the equations for the Euler method and A and b independent of time 
may be expressed as 



(26) 
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Each entry of .A is a block of the dimension of A, and each entry of x and & is a block of the dimension of x. The first 
row sets the initial value, Xq = x m . The next rows give Xj+i — (xj + Axjh) — bh, corresponding to the discretisation 
of the differential equation via the Euler method. The final rows indicate equations where Xj+i — Xj — 0. This is for 
the times where x is required to be constant. More generally, we use the Euler method as a starting method, then 
continue with a higher-order method, then for the final rows again have Xj+i — Xj — to ensure that all the final 
values are equal. Therefore, we set the blocks of A as, for N t > 2k, 



Ajj = 1, < j < k, 

i = -(1 + Ah), l<j<k, 

k+t = ail-p i Ah, k<j< N t /2, 

i = -l, N t /2<j<N t . 



3,3- 



3,3- 
3,3- 



N t /2<j<N t , 
< £ < k. 



(27) 



We will always require N t > 2k when using A, because otherwise there are not enough time steps to start the linear 
multistep method. We also set the blocks of b as 



bo = x in , 

bj =bh, 1 < j < k, 

b 3 =E k e=ofa bh > k<j<N t /2, 

bj = 0, N t /2 < j < N t . 



(28) 



We require A, b, and Xi n to be sparse, with no more than s nonzero elements in any row or column. We assume 
that the oracles are of the same form as in Ref. That is, the oracle for A is a unitary operator acting as 



A \j,£)\z) = \j,£)\z®A^) 



(29) 



Note that we use superscript to denote indexing within the block A. We also require an oracle for the sparseness, that 
locates the nonzero elements. Given a function f(j,£) that gives the row index of the ^th nonzero element in column 
j, we require a unitary oracle 



F \jJ) = \j,f(jJ))- 



(30) 
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Because A is not Hcrmitian, we require a similar oracle to give the positions of the nonzero elements in a given row. 
We also require oracles to give the values and locations of non-zero elements for b and x- m . These oracles ensure that 
the initial state corresponding to b can be prepared efficiently. Alternatively, it is also possible to consider b and x- ln 
such the efficient preparation procedure of Ref. [TtJ can be used. 

A linear system of equations can be solved using the algorithm of Ref. [lfj with complexity 0(\og(N)s 4: K 2 /cl), 
where K is the condition number of the matrix A, and ej, is the allowable error. (Note that the power of s should 
be 4, not 2 as given in Ref. ["Lot].) We use the symbol el to indicate the allowable error for the solution of the linear 
systems, which is distinct from the allowable error for the solution of the differential equation. The scaling can be 
improved if the method reported in Ref. [TlJ is used (the method is not given there, but is reported as being in a 
manuscript in preparation). The scaling reported there is 0(k 1+0 ^ log c A), which does not include the scaling in s 
or e L . 



IV. BOUNDING THE CONDITION NUMBER 



To determine the complexity in either case it is necessary to determine the value of the condition number k. To 
bound this condition number we first determine bounds on the norms of \\A\\ and ||^4 _1 ||. 

Lemma 5. The matrix A, with blocks given by Eq. (|27p . satisfies \\A\\ — 0(1) provided h — 0(1/||A||). 

Proof. To determine the upper bound on \\A\\, we express A as a sum of block-diagonal matrices, and use the triangle 
inequality. Let us define A^ to be the block diagonal matrix with all entries zero, except 



3-J 

We then have 



A%t = •AjJ-t- (3 1 ) 



-4 = $>W (32) 



1=0 

so, via the triangle inequality, 

k 

MII<EM W II- (33) 

The norm of a block-diagonal matrix is just the maximum norm of the blocks, so we find 

\\A^\\<max(l,\a k \ + \0 k \h\\A\\), 

||„4 {1} || < max(l + fc||A||, K_ x | + \p k ^\h\\A\\), 

\\A {e} \\ < \at\ + \p t \h\\A\\, K£<k. (34) 
Because we require that h = 0(l/||yl||), each of these norms is 0(1), and hence the overall norm is 0(1). □ 

Lemma 6. Suppose that the multistep method is A(a)-stable, the matrix A may be diagonalised as A — VDV^ 1 , 
and the eigenvalues of A all satisfy | arg(— Aj)| < a. Then the matrix A, with blocks given by Eq. (|27[) . satisfies 
||A _1 || = 0(N t Kv), where Ky is the condition number ofV. 

Proof. To upper bound ||A _1 ||, we use a method analogous to that used to bound the error in Ref. [l3[. As in the 
condition for Theorem [3l we assume that A may be diagonalised as 

A = VDV- 1 . (35) 

Note that A need not be Hermitian, so V need not be unitary. If we define V to be to the block matrix with V on the 
diagonal, and T> to be the matrix corresponding to A except with A replaced with D, then A = VDV^ 1 . We obtain 

ll^ll^livil-ll^ii-iiv- 1 !^^!!©- 1 !!. (36) 

To bound ||A _1 || we therefore just need to bound ||2? _1 ||. 
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The matrix V corresponds to the linear multistep solution of decoupled scalar differential equations. That is, taking 
z = V~ 1 x, the differential equation becomes N x decoupled differential equations 

zW(t) = \ j zV ] (t) + [V- 1 b]W. (37) 

The matrix T> gives decoupled linear multistep solutions of each of these differential equations. It may therefore be 
written in block-diagonal form, with each block corresponding to solution of each of these decoupled equations. The 
value of ||2? || can therefore be bounded by the maximum of the norm of the inverse of each of these blocks. 

To bound the norm of the inverse, we can take T>z = y, and determine a bound on the norm of z for a given norm 
of y. We can determine this by separately examining the uncoupled blocks in T>. For each of these blocks (labelled 
by j) we have the linear multistep equation, for m = 0, . . . , N t /2 — k, 

k 

Y,(<*i- h ^> b 2 + i=y b 2 + k- (38) 

i=0 

We also have, for the initial condition, z = j/q , and for the Euler method as the starting method with < m < re— 1, 

z® +1 -0. + h\ j )zW=y® +1 . (39) 
For the end of the simulation, we have for N t /2 < m < N t , 

^+1-4? =!&-!• (40) 
We can see that Eq. (|3"5)) is equivalent to Eq. (7.11) in the method used to bound error in Ref. fl3j . We identify 

z^+i as equivalent to e m+i in Ref. |l3j ]. and yjjjL fc as equivalent to Sf l (x m ) in Ref. [T3|. As in that method, we can 
define 

E — fz M z [j] zW) T 



As the problem is equivalent to that considered in Ref. [l3j], the result given in Eq. (7.24) of that reference holds 



A m := (y [ £ +k /(a k - hk^O, 0) T . (41) 



\E m+1 \\<M[\\Eo\\+J2\\A e \\), (42) 



where M is a constant depending only on the method. Using the definition of Af gives 



\E m+1 \\ < M \\E Q \\ + X;|y^J/|a* - h\ 3 h\ 
\ t=o ) 

^Mfll^oll+El^il/Kl]. (43) 



In the last line we have used the fact that the condition of A(a) stability means a k ■ (3k > 0, so \a% — h\j/3k\ 1 < |a*| 1 - 
For the starting method, we have used the Euler method, and the result is simpler. For the Euler method, E m and 
A TO are scalars, and are just z„ and y^+v The corresponding result is therefore, for < m < k, 

/ m— 1 \ 

Mi _l \ " u,M 



\z$\<M E i^i + Eifi+ii 

V *=o / 

m 

= M E J2\y[ j] l (44) 



Here Me is the corresponding constant for the Euler method. For the end of the simulation, we have for N t /2 < m < 

Z rn+1 ~ Z ™ ~ 2/m+lJ S0 

m 

i4 ] i<i4; /2 i+ E \yf\- ( 45 ) 

t=N t /2+\ 
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Now we can bound the norm of Eq as 



k-l 



i^oii<Ei^ ] i 



m=0 

k-l 



We can use this result to bound \zih | as, for k < m < N t /2 

\zW\<\\E m+1 _ k \\ 



<M E kY,\vf\- (46) 

[ill 



x—k 

,m 



<M\^\E4 + ^y&k\/\^\\ 

I k-l m \ 

< m M E kJ2\y l P\ + E l^'l/KI ■ (47) 



\ t=o e=k / 

For convenience we define the quantity 

M T := max.{MM E k, M/\a k \,M E , 1). (48) 

Then we find that, for all < m < N t , 



£=0 

Hence we can determine an overall upper bound on the norm of z^ as 

N t / rn 



zl$\<M T J2\y l t ] l (49) 



m=0 \£=0 

<M*N?\\yWf. (50) 



Summing over j, this then gives 

\\z\\ < M T N t \\y\\. (51) 

This result means that < MxN t , where Mt depends only on the method. This then bounds the norm of A^ 1 

as 

\\A- 1 \\=0{N t K V ). (52) 

□ 

We can now use these results to bound the condition number of A. 

Theorem 7. Suppose that the multistep method is A(a)-stable, the matrix A may be diagonalised as A = VDV^ 1 , 
the eigenvalues of A all satisfy \ arg(— A^)| < a, and h = 0(l/\\A\\). Then the matrix A, with blocks given by Eq. (|27l) , 
has condition number k = O(NtKv), where Ky is the condition number of V . 

Proof. The condition number is given by the formula 



_/ \\M\\ 

K ^JxTJ \-T~\\AH 



max-p^q- I = \\A\\ ■ (53) 



The conditions of this theorem ensure that the conditions of Lemmas [5] and [5] hold. Therefore we can use the bounds 
on ||^4|| and ||^4 _1 || from those lemmas to obtain k — 0(N t Kv)- D 

[7*1 

This result for the condition number can be explained in a more intuitive way. Each value of y. is equivalent to an 
excitation of the differential equation at a single time. Therefore z is close to the solution of the differential equation 
with each of those individual excitations. An excitation can not cause a growing solution, because of the stability 
condition. This means that the worst that an excitation can do is cause a solution that is displaced by a proportional 
amount for the remaining time. Therefore the norm of z can not be more than a factor of Nt times the norm of y. 
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V. ALGORITHM FOR SOLVING LINEAR SYSTEMS 

We can use the bound on the condition number to estimate the complexity of the quantum algorithm. We will first 
explain the scaling in a simple way, then give a more rigorous result. Using Eq. (|24[) for the number of time steps, we 
have (ignoring dependence on many of the quantities) 



(\\A\\At) 



l+i/p 



K = \^A )■ ( 54 ) 

Using this expression in the result for the complexity of solving linear systems from Ref. fiol ] gives 

d{\og{N) S \\\A\\M) 2+2 lv/{e 2 IPe L )). (55) 
If the technique of Ref. can be used, then the scaling can be improved to 

0(log c (AO(PI|Ai) 1+1/ 7e 1/p ). (56) 

There is a question of whether the scaling in Ref. [Io| can be improved because we consider a specific application 
of the solution of linear systems. The scaling obtained in Ref. [l(| is based on a worst-case scenario, that scales 
as A max ||x||. It is easily shown that ||x|| scales as \/Nt, provided the magnitude of the solution of the differential 
equation does not vary greatly in time. In contrast, the magnitude of b is given by 



l&H 2 = \\x iD \\ 2 + (k- l)\\b\\ 2 h 2 + (N t /2 -k + l)\\b\\ 2 h 2 (X>J 

< || 2 ; in || 2 + (fc-l)||&|| 2 / l 2 + /iAt||&|| 2 • ( 57 ) 



si=0 



In the case that we use the Euler method, then h oc 1/Ai, so ||6|| = 0(Vl In that case it is possible to solve the 
linear systems more efficiently than the worst-case scaling given by Ref. [1(|. However, here we are concerned with 
using higher-order linear multistep methods. For these methods the scaling of hAt is close to N t . In that case ||6|| 
has similar scaling to A max ||x||, so there is not a significant advantage to using this approach to improve on the result 
in Therefore we will just use the result stated in Ref. [Icj . 

Before obtaining the overall scaling, another factor that needs to be considered is the scaling for creating the state 
encoding b. This is because the algorithm of Ref. [10( uses an amplitude amplification approach, which requires the 
preparation of this state at each step. Therefore, the overall complexity is multiplied by the complexity of performing 
this state preparation. We assume that we have oracles that give the elements of x- m and b in the computational basis, 
and that these vectors are s-sparse. We find the following result. 



Lemma 8. The state encoding b, 



(58) 



can be prepared using 0(y/s + log(A^)) calls to the oracles for x- m and b, provided the normalisations of x ln and b are 
known. 

Proof. For this preparation we can assume that the normalisations of x- m and b are known. That is because this 
normalisation can be determined with complexity O(s), and need only be determined once. Because the overall 
complexity of the simulation is greater than this, it can just be assumed that the normalisation is known. 

A state of dimension s can be prepared with complexity 0(yfs) using the method of Ref. (18(. As discussed above, 
we assume an oracle of the form used in Ref. |8| for the sparseness. This means that it only requires one oracle call 
to prepare an s-sparse state from a dimension s state with the same coefficients Q. Therefore the complexity of 
preparing an s-sparse state is also 0(yfs). 

Let us encode the state in three registers (which may themselves be composed of multiple qubits). The first is one 
qubit encoding |0) for the times ti,... ,tjy t /2, and |1) for the times fo,iiv t /2+ij The second register provides 

the remainder of the encoding of the time. The third register is of dimension N x , and encodes b or x- ln . 
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By performing a rotation on the first qubit, based on the normalisations of cc; n and b, we obtain the correct relative 
weighting of the two time ranges. Then, conditional on the qubit being in the state |0), we can prepare a superposition 
of times ti, . . . ,tjj t /2 m the second register, as well as a state encoding b in the third register. Conditional on the 
qubit being in the state |1), we can prepare the time to in the second register, and x- ln in the third register. The 
complexity of these controlled state preparations will be 0(y / s). The complexity of preparing the superposition of 
times can be made 0(\og(N t )) simply by choosing N t to be a power of two (which does not change the scaling). □ 

We now translate the complexity of solving linear systems into the complexity of obtaining a state corresponding 
to the solution of the differential equation. 

Theorem 9. Suppose that the multistep method is order p and A(a)- stable, the matrix A may be diagonalised as 
A = VDV" 1 , the eigenvalues of A all satisfy | arg(— Aj)| < a, and 

max \\x{t)\\=O{\\x(t + At)\\), (59) 

*6[*o,*o+At] 

e = o(||x in ||). (60) 

Then a state encoding the solution of the differential equation at time to + At may be obtained to within trace distance 
e using 

6 (log(N x y/ 2 (\\A\\At)V\l 3/ \\\x^ + ||&||/||A||) 5 / 4 ||*(t + At||/e 9 / 4 ) (61) 
calls to the oracles for A, b, and X{ n . 

Proof. There are two main issues that we need to take into account in determining the complexity of obtaining the 
solution of the differential equation. The first is the probability of obtaining the correct time, and the second is the 
relation of the allowable error for solving the differential equation to the allowable error for the solution of the linear 
systems. 

To obtain the correct final state, one would need to measure the time in the ancilla register, and obtain a value in 
the range to + At to to + 2At. This probability can, in principle, be small or even zero. However, the conditions of 
the theorem ensure that it is not. The probability of success is given by 



E 



j=N t /2 \\ X 3 



pti ™ = — m — ■ (62) 



The normalisation of the state is given by 



N t 

£ini 2 

.7=0 

Ei 

3=0 



x(t + hj)\\ + 0(e)} 2 

= 0[N t max ||x(t)|| 2 
V te[t ,*o+At] 

= O (N t \\x(t + At)\\ 2 ) . (63) 

Here we have bounded the error in the state at all times by e. This is because we choose parameters that bound the 
error at time to + At by e. The bound on the error increases monotonically with time, so the error at earlier times 
will also be bounded by e. We have also used the conditions (|59| and (|60|) to obtain that e = o(||x(to + At)||). Using 
this result for (tp\ ip) we obtain 



Ptimc ^ 



'Y^N t /2[Mto + At)\\+0(e)] 



N t \\x(t + At)|| 2 J 

= m- (64) 

Therefore the probability of obtaining the correct time does not change the scaling of the complexity. 

The other main issue that needs to be addressed to obtain the overall scaling is the scaling in the error. The €l 
used in the expression for the solution of the linear systems is the error in the state encoding all times, not the error in 
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the estimate of a; (to + At) obtained (for which we use e). To determine the relationship between e and cl, we use the 
bound on the normalisation of the state. In solving for the state we obtain the state coefficients approximating 



x(t Q + AfJ/VW)- (65) 

Using the bound on the normalisation in Eq. (|63|) . if the error in the coefficients is no more than el, the error in 
the solution of the differential equation will be 0(eLy/Nf\\x(to + Af||). Therefore, the error in the solution of the 
differential equation will be bounded by e if we take 

e L = e(e/[ v / A 7 Mt + At)||]) . (66) 

Using this value of cl in the scaling from Ref. [ioj . together with k — 0(N t Kv) from Theorem^ and the complexity 
of state preparation from Lemma |8l the number of oracle calls is 

6(log(N x ) S 9 / 2 N^\ 2 v \\x(t + At\\/e). (67) 

Note that we can omit \og(N t ) from Lemma[8j because the O notation omits logarithmic factors. For the same reason, 
we have replaced N = N x N t with N x . We use a value of N t that is sufficient to ensure that the error is no greater 
than e, and that h = 0(1/\\A\\), which is a condition needed to use Theorem [7] If we take 



n, = e I iiaiia^^MIM + h/pii> + mm ^„ ± w/Mlj y m 

then using LemmaSl the error will be bounded by e. In addition, because e = o(||xi n ||), we obtain N t = f2(||A||At), 
which ensures that h = 0(1/\\A\\). 
We can simplify the result by taking 



N t = Q ^(pnA^i-n/P^dl^H + 11^11/11^11) j . (69 ) 

Then the overall scaling of the number of black-box calls is 

6 (log(Ag s 9 / 2 (||A|| At)( 5 / 2 )( 1+1 /f)4 3/4 (||x in || + \\b\\/\\A\\) 5 ^\\x(t + At\\/e 9 / 4 ) . (70) 

Because we are using the O notation, which omits sublinear terms, we can omit l/p in giving the result. □ 

This result is somewhat conservative, because we have included the term due to the error in starting the linear 
multistep method. If we assume that this error is negligible, then we obtain 

6 (log^j^ciiAiiAtj^/^a+i/pj^+Vp^^ii + ii^i/pid^ii^ + At||/ e 1+5 / 2 p) . (7i) 

This has the same scaling in ||^4||A£, but improved scaling in other quantities. 

VI. CONCLUSIONS 

A quantum computer may be used to solve sparse systems of linear differential equations, provided the result may 
be encoded in a quantum state, rather than given explicitly. By encoding the differential equation as a linear system, 
and using the algorithm of Ref. 10] for solving linear systems, the complexity is (including only scaling in ||j4|| and 
At), 

6 ({\\A\\Mfl 2 \ . (72) 



This improves upon previous results for nonlinear differential equations, which were strongly exponential in At 0]. 
This algorithm has an enormous range of possible applications, because large systems of differential equations are 
ubiquitous in science and engineering. In particular they arise from the discretisation of partial differential equations. 
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These results are for constant coefficients, because that enables an analytic error analysis. This approach can also 
be used to solve linear differential equations with time-dependent coefficients, though the error analysis will be more 
difficult. 

An interesting question is how the complexity will scale if the method referred to in Ref. [Ill ] can be used. That 
reference refers to work in preparation providing an improved scaling for solving linear systems. The scaling is close 
to linear in the condition number, which would improve the scaling in ||_4||Ai for the solution of linear differential 
equations. 

Another interesting direction for future research is the problem of nonlinear differential equations. It is likely that 
the exponential scaling obtained in Ref. Q is fundamental, because quantum mechanics is linear. However, it may 
potentially be possible to improve the constant of the exponential scaling. 
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